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ABSTRACT 

We model a spiral galaxy by a thin axially symmetric disk that includes both 
^ I visible and dark matter. The surface mass density of the disk is calculated directly 

Q\^ I from the rotational velocity curve without extra assumptions. We simplify the 

O ! to some radius, Vmax, we extrapolate them by attaching a Keplerian tail. The 



O 



^ 



a 



numerical procedure and the extrapolation are tested with a known toy mass 

density and shown to reconstruct it with a good precision if Vmax includes a 

Qh| sufficient part of the velocity curve. Mass density curves are calculated for Milky 

O I Way and NGC 3198. We vary the extent of the flat part of the velocity curves 

^ I from 30 kpc to 200 kpc and show that does not affect appreciably the calculated 

mass density inside Vmax = 30 kpc . The reconstructed masses for Milky Way are 

15 X 10^° Mq inside the visible disk and 23 x 10^° Mq inside 30 kpc. For NGC 

/\ ' 3198, the reconstructed mass inside the visible disk is 6.5 x lO^'' Mq and 11 x 

c3 I 10^° Mq inside 30 kpc. The total galactic masses are roughly proportional to the 

extent of the fiat part of the velocity curves which is currently unknown. The high 

light-to-mass ratios obtained for the visible disks of the galaxies - 11 Mq/Lvq for 

Milky way and 9.3 Mq/Lyq for NGC 3198 - suggest presence of dark matter. The 

method is also applied to NGC 3031 - a spiral galaxy with a declining velocity 

curve in which case it is able to reconstruct both the mass density curve and the 

total mass (14 x 10^° Mq). 

Subject headings: dark matter — galaxies: halos — galaxies: kinematics and 
dynamics 
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INTRODUCTION 



The flat rotational velocity curves of spiral galaxies remain a central problem in as- 
trophysics. Exotic attempts for explanation were made in terms of Modified Newtonian 
Dynamics (Milgrom 1983) or even incorporating general relativity effects (Cooperstock and 
Tieu 2005). The most popular approach remains a three component galaxy model consisting 
of a bulge, a thin disk, and a spherical dark matter halo (Kent 1987). Usually, the mass- 
to-light ratios of the bulge and the disk are taken constant and a the halo mass density is 
assumed to vary as 1/[1 + {r/aY] since that produces a fiat velocity curve by itself. The sev- 
eral parameters of the model are then fit to the observed velocity profile. Unfortunately, such 
fits do not produce unique values of the parameters and frequently some of them have to be 
chosen before the fit (Binney and Merrifield 1998). A model with fewer a priori assumptions 
is desirable. 

In this paper, a spiral galaxy is represented by an axially symmetric thin disk of surface 
mass density a{r) which includes all types of matter. It is very natural to assume that the 
dark matter component (if any) resides in the galaxy disk just like the visible matter since 
both are supposed to have the same gravitational properties. The advantage of that model 
is that a{r) can be calculated directly from the velocity curve v{r) without invoking extra 
assumptions. The necessary formulas were developed by Toomre (1963). Since the velocity 
curves are known out to a radius r^ax) they have to be extrapolated by a Keplerian decline 
(Nordsieck 1973). The thin disk model is applicable in cases where the galaxy radius is much 
bigger than its thickness and the size of the bulge. It ignores deviations from axial symmetry 
due to the spiral arms or possible bars. Newtonian gravity law is used throughout assuming 
the general relativity effects are negligible which is the prevalent opinion in the field. The 
model does not address the question why the mass is distributed in the way calculated. 

In the current paper, we apply the model in a way that in our opinion is simpler than 
the standard approach. Using a toy mass density, we test numerically how the accuracy of 
the reconstructed mass density depends on the point, r^axi after which the velocity profile 
is extrapolated. At the end, the density curves of several galaxies are calculated using 
reasonable guesses of the extent of the fiat part of the velocity curves. Luckily, the density 
curves and the masses inside the current experimental Tmax are not very sensitive to that 
extent. On the other hand, only a range of the total masses of the galaxies can be given, 
since they are roughly proportional to the unknown extent of the velocity curve flatness. 
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2. RELATION BETWEEN MASS DENSITY AND VELOCITY CURVE 

FOR A THIN DISK 

Outside the rotating disk, the gravitational potential satisfies the Laplace equation 

V2$ = (1) 

whose general solution in cylindrical coordinates for 2; > is (Jackson 1998, p. 118) 

<l>(r,0,z) = y^/ dk e~''''Jm{kr)[Am{k) COS m(f) + Bm{k) sin m(f)] (2) 

where Jm are the usual Bessel functions of the first kind and Am and Bm are expansion 
coefficients. The solution for 2; < is just a mirror image of the above. An axially symmetric 
solution is given by setting m = 0: 

$(r,2)=/ dke-^'Jo{kr)A{k) (3) 

Jo 

where A{k) = Ao(^) denote the Bessel spectrum of the potential. Assuming a finite mass 
density, the gravitational potential must be continuous on the disk plane, z = 0, and is given 
by 

$(r, z = 0)= dk Jo{kr) A{k). (4) 

Jo 

This equation can be inverted 



A{k) = k drr Jo{kr) <l>(r, z = 0) (5) 

Jo 

using the orthogonality relation (Jackson 1998, p. 118) 

/ dr rJm{kr)Jm{k'r) = -6{k - k'). (6) 

Jo k 

The rotational velocity profile, f (r), is connected to the potential by 
t;2(r) 9<l'(r,0) 



dr 



dk k Ji{kr)A{k) (7) 



where we used dJo{x)/dx = —Ji{x). This can be inverted using equation (6) with m=l 
giving 



00 

2/ 



A{k) = - / dr Jiikr) v\r). 
'0 
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On the other hand, Gauss theorem connects the surface mass density, o"(r), to the 
potential by 



2 d^(r,z) 
a r) - 



AivG dz 



= 7^ rdk{-k)Mkr)A{k), (9) 



where G is the gravitational constant. The Bessel functions Jo{kr) form a complete set and 
the mass density can be expanded analogously to the potential 

a{r) = dk Joikr) a{k) (10) 

where <j{k) is the Bessel spectrum of the density. Comparing the last two equations gives 
the relation between the Bessel spectra of the potential and the surface mass density 

.W^-H^. (11) 

The way to compute the mass density for a given velocity profile, f (r), is now clear: 

1. Use equation (8) to compute A{k), the Bessel spectrum of the gravitational potential. 

2. Obtain the Bessel spectrum of the mass density, (y{k) from equation (11). 

3. Reconstruct the surface mass density, o"(r), from equation (10). 



3. EXTRAPOLATING THE VELOCITY CURVE 

The velocity curves of galaxies are known from HI spectroscopy out to some radius, 
"^max) of about 20 — 40 kpc. They have to be extrapolated to r = oo before computing 
equation (8). 

The most elementary approach used in Takamiya and Sofue (2000) would be to perform 
the integration in equation (8) up to Vmax ■ That is equivalent to setting the rotational speed 
to zero for bigger radii. Unfortunately, as the calculated a{r) has shown, that produces a 
ring of negative mass density around r^ax which is necessary to cancel out the field of the 
positive mass so that v{r) can drop to zero at that radius. Negative mass density is hard to 
justify physically. Another possibility is to fit the velocity profile and let the continuation 
of the fitting function determine the behavior for r > Vmax- That is not acceptable either 
because the fit continuation usually does not have the necessary Keplerian decline at bigger 
radii. 
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If we assume that Vmax is close to the Keplerian tail of the velocity curve, we can estimate 
the total mass of the disk, Me, using the point-mass formula 

That estimate is good if all or a big part of the mass is inside Vmax- In this paper, we 
extrapolate the velocity curve by attaching a Keplerian tail to it: 



.„,,) . ^^. , > ,_. (13) 

a similar approach was used by Nordsieck (1973). 

Numerical integration of Bessel functions over a large range like in equation (8) can give 
severe cancellation errors. That is why all authors (Toomre 1963; Nordsieck 1973; Takamiya 
and Sofue 2000) avoided the simple Bessel spectrum formulation presented in the previous 
section and would rather integrate over the velocity derivative dv^{r)/dr and use elliptic 
integrals. That has the inconvenience that the elliptic integrals are divergent at a point and 
the velocity derivative has to be computed numerically from the velocity curve which is a 
possible source of errors. The success of our computations shows all these complications 
are not necessary since we take the integral of equation (8) numerically up to r = fmax and 
analytically after that: 

/~ * Mkr) vHr) = GM, f *■ ^^M . GM^ P " *^ ^1^' {§■ ^J- -(^)^)1 • * ^ » 

•Jrmax J rmax ^ L U , /C — U 

(14) 
Here -F[- ■ ■] denotes the generalized hypergeometric function for which good numerical rou- 
tines exist. It is defined by the expansion 

oo T{a+n) ^ 

n=0 r(bi) T{b2) 

where r(a;) denotes the usual gamma function. 



4. CALCULATION PROCEDURE 

The velocity curve is a smooth fit to the experimental points out to r^nax- It is extrap- 
olated according to equations (12) and (13). The total mass of the computed a{r) will be 
Me since it is already encoded in the extrapolated Keplerian tail of v{r). The gravitational 
potential Bessel spectrum, A(k), of equation (8) is sampled for k = 0...kmax with a step of 
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AA;. The corresponding mass density spectrum, (T{k), is computed from equation (11) and 
interpolated with sphnes. The values for k^ax and AA; are chosen so that the computed 
spectrum is without sudden kinks and stabilizes to zero before kmax- The spline-fitted a{k) 
is integrated numerically according to equation (10) returning a{r) at the points of interest. 



TESTING THE METHOD WITH A TOY MASS DENSITY 



In this section we work in arbitrary units and the gravitational constant G = 1. Our 



toy surface mass density shown on Fig. 1(a) is 



a(r) 



' 5.5 + 0.5 cos(27rr) ,0<r<2 
^ ,2 <r <4 

,4<r < 6 
,6<r < 8 
,r > 8 



l+(r-2)2 

1.2 

4.8- 0.6 r 
U 



(16) 



The chosen a{r) is continuous since the gravitational potential and the rotational speed 
would diverge at points of discontinuity. That is directly related to the assumed zero thick- 
ness of the disk. The form of a{r) does not have a physical significance. Its solely purpose is 
to test the ability of the method to reconstruct it from an incomplete velocity profile. The 
constants in o"(r) are chosen so that the radial derivative of the produced potential is always 
positive i.e. the gravitational force is towards the center not away from it. The piece for 
< r < 2 will test the method's ability to reconstruct small details in o-(r), the decreasing 
piece for 4 < r < 6 is a form that produces a flat velocity curve by itself, the piece for 
4 < r < 6 will test how the method reconstructs a planar halo, and the piece for 6 < r < 8 



original mass density 





10 12 



Fig. 1. — (a) original (test) mass density (b) integration variables in the disk plane. 
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is just a continuous drop to zero. 

The total mass is 

M 



a{r)2Txrdr = 308.6. 



;i7) 



The quantity relevant to the amount of mass within certain radii is ra{r) not a{r) alone. 
That amplifies the mass density at bigger radii and the 4 < r < 8 part of the mass density 
actually contains 40% of the total mass although it looks negligible on the a{r) plot. 

Before testing the method, it is necessary to compute the potential and the real velocity 
curve from the given toy mass density. It is possible to use the derived Bessel spectrum 
formulas in reverse mode but we chose an independent method - direct numerical integration. 
The gravitational potential at radius r on the disk is (see Fig. 1(b)) 

f -Gdm /"^^ , „ ^.., , 

$(r,z = 0)=/ = -Gl / . ' ' " . fl^ 



P J4,=o Jr=o ^R^ + r2 - 2Rrcos(f)' ' ' 

This form is problematic for numerical integration at the points where the denominator of 
the integrand goes to zero. That zero is canceled by the numerator but most numerical 
routines are likely to have severe cancellation errors. The problem is eliminated if dm is 
expressed in terms of coordinates, {p,0), with respect to the point where the potential is 
being computed (see Fig. 1(b)) since that cancels out the denominator: 



^{r,z = 0) 



-Gdm 



-G 



P 



-G 



61=0 J p=0 



P 



a{R)dpde, (19) 



e=o Jp=o 



where R = ^Jr"^ + p^ — 2 r p cosQ. 

gravitational potential 



velocity curve 





Fig. 2. — toy mass density: (a) gravitational potential (b) velocity curve. 

Using equation (19), the potential was calculated numerically for radii from to 20 
with a step of Ar = 0.1 and interpolated with splines (shown on Fig 2(a)). The derivative 



of the spline-fitted potential was computed numerically and then the velocity curve (shown 
on Fig. 2(b)) was calculated according to 



f (r) 



d^{r,z = 0) 



dr 



(20) 



In the following, we test how well the method will reconstruct the toy mass density if 
part of the velocity curve is not known, hence extrapolated. 

In the first test we assume that the velocity curve is known out to r^ax = 10. The 
total mass estimate obtained from equation (12) is Me = 368.2 which overshoots the real 
mass M = 308.6 The velocity curve was extrapolated after Vmax (see Fig. 3(a)) by formula 
(13) and fed to the calculation procedure outlined in the previous section. The resulting 
reconstructed mass density is shown on Fig. 3(b). It is so close to the original one that 
they are hard to distinguish. The original and the reconstructed integrated masses for 
^ < f^max are, respectively, 308.6 and 319.1. The higher reconstructed mass is caused by the 
overshooting of the total mass estimate Mg. 



original and extrapolated velocity curve 



original and reconstructed mass density 





Fig. 3. — (a) original (dashed) and extrapolated (solid) velocity curve for r^ax = 10 (b) 
original (dashed) and reconstructed (solid) mass density. 

In the second test we assume that the velocity curve is known out to Vmax = 4. That 
leaves 40% of the original mass outside Vmax- The total mass estimate is Mg = 226.0 which is 
below the real mass M = 308.6. The extrapolated velocity curve (see Fig. 4(a)) has lost one of 
the 'humps' and is below the original one for r > Vmax- Due to that we expect a reconstructed 
mass density below the original one and a poor reconstruction of the mass density shape for 
r > Tmax- The reconstructed mass density shown on Fig. 4(b) indeed preserves the shape of 
the original one although it is slightly below it for r < Vmax- The reconstructed mass inside 
f^max is 158.5 while the original mass inside Vmax is 182.9. The discrepancy again is explained 
by the low mass estimate M^. 
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original and extrapolated velocity curves 



original and reconstructed mass density 
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Fig. 4. — (a) original (dashed) and extrapolated (solid) velocity curve for r^ 
original (dashed) and reconstructed (solid) mass density. 



4 (b) 



In the third test we assume that the velocity curve is known out to Vmax = 2. That 
leaves 61% of the original mass outside Vmax- The total mass estimate is Mg = 50.6 which is 
far below the real mass M = 308.6. The extrapolated velocity curve (see Fig. 5(a)) has lost a 
lot of details and height for r > Vmax- Due to that we expect the reconstructed mass density 
to be severely below the original one and poor reconstruction of the mass density shape for 
r > Tmax- The reconstructed mass density shown on Fig. 5(b) indeed preserves somewhat the 
shape of the original one although is very much below it for r < Vmax- The shape is totaly 
washed out for r > Vmax- The reconstructed mass inside Vmax is 27.4 while the original mass 
inside r^ax is 69.1. 



original and extrapolated velocity curve 



original and reconstructed mass density 
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Fig. 5. — (a) original (dashed) and extrapolated (solid) velocity curve for r^ 
original (dashed) and reconstructed (solid) mass density. 



2 (b) 



The tests show that the total mass estimate M^ is close to the real mass when Vmax is 
on the Keplerian tail of the velocity curve (test 1) or at least near it (test 2) i.e. when a big 
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part of the total mass is inside Vmax- The height of the reconstructed mass density depends 
on Me but is not extremely sensitive to it - Mg is off by 20% in test 1, yet the reconstructed 
distribution is indistinguishable from the original. The shape of the reconstructed mass 
density inside Vmax is fairly independent of the goodness of the mass estimate. The shape 
outside Tmax is washed out due to the velocity curve extrapolation (tests 2 and 3). 



6. APPLICATION TO SPIRAL GALAXIES 



The velocity curve fits of several galaxies are taken from Cooperstock and Tieu (2005). 

The velocity curve of Milky Way is known from HI observations out to Vmax = 30 kpc. 
The total mass estimate of equation (12) is proportional to the extent of the flat part of 
the velocity curve, which is currently unknown. We can speculate that the Keplerian tail 
starts right after 30 kpc, which gives a total mass estimate of Mg = 30 x 10^*^ Mq. The 
corresponding extrapolated velocity curve and the reconstructed mass density are shown on 
Fig. 6. The visible disk has a radius of 15 kpc. The reconstructed mass is 15 x 10^*^ Mq (half 
of the total mass) inside the visible disk and 23 x 10^° Mq inside Vmax- This is in contrast to 
the widespread view that the mass must be concentrated predominantly outside the visible 
disk. The conventional value of the Milky Way's absolute magnitude is -20.5 which gives a 
luminosity of Ly = 1.36 x lO^^LyQ. The mass-to-light ratio is M/Ly = 11 Mq/ LyQ for the 
visible disk and 22 Mq/ LyQ for the whole galaxy. 



. q5 extrapolated velocity curve 



reconstructed surface mass density 





15 20 25 
r, Kpc 



Fig. 6. — Milky Way: (a) velocity curve extrapolated after rmax = 30 kpc (b) reconstructed 
mass density for r^nax = 30 kpc (solid) and r^ax = 200 kpc (dashed). 



To evaluate the effect of the unknown extent of the fiat part of the velocity curve, we 
continued the experimental curve with a flat line till r^ax = 200 kpc and extrapolated with 



a Keplerian tail after that. That gives a total mass estimate Mg = 197 x 10 Mq. The 
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resulting reconstructed mass density is the dashed curve on Fig. 6(b). The reconstructed 
mass inside the visible disk is 16 x 10^° Mq and inside 30 kpc is 23 x 10^*^ Mq. It is remarkable, 
that although the total mass estimate increased more than 6 times, the masses and the mass 
density inside the experimentally known part of the velocity curve (r < 30 kpc) remained 
fairly unchanged. 

Another well studied and measured galaxy is NGC 3198 whose velocity curve is known 
out to r^ax = 30 kpc. Assuming the Keplerian decay starts right after r^ax, gives Mg = 
15 X 10^*^ Mq. The extrapolated velocity profile and the resulting mass density are shown 
on Fig. 7. The visible disk has a radius of 14 kpc. The reconstructed mass inside the visible 
disk is 6.5 x 10^° M© (43% of the total mass) and 11 x 10^° Mq inside Vmax ■ The V-band 
luminosity of NGC 3198 is Ly = 0.7 x lO^^Ly© (van Albada et al. 1985). The mass-to-light 
ratio is M/Ly = 9.3 Mq/Lvq for the visible disk and 21 Mq/Lvq for the whole galaxy. 

Continuing the fiat part of the velocity curve with a straight line till Vmax = 200 kpc 
and extrapolating with a Keplerian tail after that gives Me = 100 x 10^° Mq. The resulting 
mass density is the dashed line on Fig. 7(b). The corresponding mass inside the visible 
disk is 7.2 x 10^° Mq and 15 x 10^° Mq inside 30 kpc. Again, the mass distribution inside 
f^max = 30 kpc is relatively insensitive to the extent of the velocity curve fiat part. 
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Fig. 7. — galaxy NGC 3198: (a) velocity curve extrapolated after r^ax = 30 kpc (b) recon- 
structed mass density for r^ax = 30 kpc (solid) and r^ax = 200 kpc (dashed). 



Last, we consider the galaxy NGC 3031 whose velocity curve is known out to rmax = 
21 kpc. It is an excellent candidate for extrapolation since it seems to show a Keplerian 
decline already. The extrapolated part of the NGC 3031 velocity curve continues very nat- 
urally the slope of the experimental curve (see Fig. 8(a)). We can expect an excellent 
reconstruction of the mass density curve possibly with some features outside r^ax washed 
out. The resulting surface mass density is shown on Fig. 8(b). The total mass estimate is 
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Me = 14 X 10^° Mq. The reconstructed mass inside r^nax is 11 x 10^° Mq (79% of the total 
mass). 
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Fig. 8. — galaxy NGC 3031: (a) velocity curve extrapolated after r^ 
structed mass density. 
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DISCUSSION 



For a spherically symmetric mass distribution, the mass density inside some radius can 
be calculated without knowing the velocity curve outside that radius. In the axial symmetry 
of the disk model, that is generally not true. Nevertheless, the reconstructed mass densities 
inside r = 30 kpc were shown in the previous section to be rather insensitive to the extent of 
the flat part of the velocity curve which was varied from r^^ax = 30 kpc to 200 kpc. We can 
regard the obtained mass densities and the total masses inside r^ax = 30 kpc as fairly good 
estimates of the real ones unless the model is not applicable or the velocity curves start to 
increase after r^ax- 

The high mass-to-light ratios obtained, even when the velocity curves were assumed to 
decay after the current experimental Vmaxi are suggestive of the presence of dark matter. 
What fraction of the matter is dark or nonbarionic remains to be seen. 

The estimates of the total galactic masses are roughly proportional to the extent of the 
flat parts of the velocity curves, which are currently unknown. While the mass inside the 
visible disk remains relatively constant, the mass outside it increases with increasing the 
total galactic mass. 

If applicable, the fiat disk model with velocity extrapolation, is expected to give fairly 
accurate estimates both of the total mass and the surface mass density function for galaxies 
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with already decaying velocity curves like NGC 3031. 
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